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We present detailed simulations of a generalization of the Domany-Kinzel model to 2+1 dimen- 
sions. It has two control parameters p and q which describe the probabilities Pk of a site to be 
wetted, if exactly k of its "upstream" neighbours are already wetted. If Pk depends only weakly 
on k, the active/adsorbed phase transition is in the directed percolation (DP) universality class. If, 
however, Pk increases fast with k so that the formation of inactive holes surrounded by active sites 
is suppressed, the transition is first order. These two transition lines meet at a tricritical point. 
This point should be in the same universality class as a tricritical transition in the contact process 
studied recently by Liibeck. Critical exponents for it have been calculated previously by means 
of the field theoretic epsilon-expansion (e = 3 — d, with d — 2 in the present case). Rather poor 
agreement is found with either. 



I. INTRODUCTION 

The critical behaviour of directed percolation (DP) has 
been studied since more than 30 years, if we also count 
the work on 'reggeon field theory' |l| which was only later 
recognized as a field theory for a stochastic process which 
is in the DP universality class [2|. In spite of this long 
history, several of its aspects have barely received any 
attention yet. In particular this is true for its tricritical 
version. 

DP is e.g. realized by a reaction-diffusion system 



A -> 2 A (rate kx), 
A -> (rate k 2 ), 
2 A -> A (rate k 3 ). 



(1) 



Its mean field theory is just the rate equation for the 
number n(t) of A particles, h = an — bn 2 with a = k\ — k 2 
and b = k 3 . The transition is in this description at a = 
0. The first mentioning of tricritical DP (TDP) seems 
to have been made in this context by Q, where it was 
pointed out that the system composed of the reactions 



2 A -> 3 A (rate k[), 
2 A -> (rate k' 2 ), 
3 A ->• 2 A (rate k' 3 ) 



has the rate equation 



= bn 2 - , 



(2) 



(3) 



with 6 = k[ — 2k' 2 ,c — k' 3 . This is just the mean field 
equation for a system with tricritical point at 



b = 0, c> 0. 



(4) 



But, as observed by Janssen Q and Ohtsuki & Keyes 
Eqs.@ and (@J can be realized also differently, e.g. by 
combining both reaction-diffusion systems and choosing 
ki-k 2 = 0, k[ - 2k' 2 -k 3 = 0,k' 3 > 0. 

Both versions differ however beyond the mean field 
level, if fluctuations are taken into account. While the 



combination of Eqs.QJ and (J2J should be described by 
a Langevin equation n(x, t) — an(x, t) + 6n 2 (x, t) — 
cn 3 (x, t) + £(x, t) with the "semi-multiplicative" noise, 
(£(x,*)£(x',t')) oc n(x,t)S(x ~ x')S(t ~ t'), well known 
from reggeon field theory, the noise appropriate for Eq. J2J) 
alone should be particle number conserving when n — > 0, 
i.e. @ (£(x,i)£(x',t')) oc n 2 (x,i)(5(x - x')S(t - t') + 
const n(x, t)V 2 S(x - x')5(t - t'). 

The reaction scheme Eq. has become infamous dur- 
ing the last years as "pair contact process with diffu- 
sion" (PCPD) pfl, and has received very much atten- 
tion. It is still basically unsolved 0. In contrast, there 
was no activity on the "true" tricritical DP with semi- 
multiplicative noise, after Ohtsuki & Keyes had worked 
out the lowest order results of the e— expansion. This 
has changed only very recently with field theory work by 
Janssen 191 and with extensive and careful simulations by 
Liibeck [lOj . 

The present work was mainly started because only 
the stationary behaviour was studied in [ljj, and we 
wanted to obtain also dynamical (tri-) critical exponents. 
It seems that the dynamical behaviour of TDP was never 
studied numerically before. Indeed, for ordinary DP it is 
also easier to obtain static critical exponents from dy- 
namical simulations than from stationary behaviour. Al- 
though this might be different for TDP, it seems worth 
to check it. But we made also a few other changes. In 
particular we study a model with parallel updating in 
discrete time, while a model with continuous time (the 
tricritical 'contact process' 01 instead of TDP proper) 
was studied in [Io|. These two models should be in the 
same universality class, but a check would be welcome 
- and simulating discrete time processes is often faster. 
Finally, while the exponents observed in |10| were rather 
close to the mean field values, some of them deviated from 
mean field theory in the opposite direction of that pre- 
dicted by the e— expansion. Simulating a different model 
in the same universality class with emphasis on different 
aspects might clarify whether this indicates a failure (or 
slow convergence) of the e-expansion, or hints at prob- 
lems with the simulations. 

In this paper, we study a fully discrete model which 
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is actually a generalization of the well known Domany- 
Kinzcl (DK) model 12]) to 2+1 dimensions. This gener- 
alization to higher dimensions is necessary, because there 
is no non-trivial tricritical point in 1+1 dimensions (the 
cross over to the first order transition, called "compact 
DP" [lj| in this case, happens at the boundary of the 
DK phase diagram, just as the phase transition in the 1-d 
Ising model occurs at T — 14]. At the tricritical point 
we obtain rather precise estimates of the (tri-)critical ex- 
ponents, except for the cross-over exponent which is af- 
fected by large corrections to scaling. Beyond this point, 
when the active/ adsorbed transition is first order, we 
find that clusters starting with single sites survive with 
a stretched exponential probability. This is similar to 
the decay of clusters in the Grassberger-Chate-Rousseau 
(GCR) model with re-infection easier than first infection 
|l5l Il6| , and has a similar qualitative reason |Tt| . The 
boundary between active and adsorbed regions behaves, 
near the first order transition, like a generic fluctuating 
non-equilibrium (KPZ |18() surface and shows the same 
scaling laws. 

In the next section we will define the model and de- 
scribe some special features of the simulations. Predic- 
tions to compare our simulations with are reviewed in 
Sec. 3. Results will be presented in Sec. 4, while the paper 
end with a discussion in Sec. 5. 



II. THE GENERALIZED DOMANY-KINZEL 
MODEL 

The standard DK model lives on a square lattice. Usu- 
ally, this lattice is drawn as tilted by 45°, and a site can be 
wetted by its two lower neighbours. For a more easy later 
generalization to 2+1 dimensions, we use a non-tilted lat- 
tice, so that site (x,t + 1) can be wetted by sites (x — l,t) 
and {x + 1, t). The probability to be wetted is P\ = y, if 
one of these sites was wetted, and P% = z, if both were 
wetted. Here y and z are real numbers between and 
1. If y = z, this corresponds to site DP, while bond DP 
corresponds to z — y(2 — y). If y < 1/2, any finite cluster 
of active (=wetted) sites dies with probability one. For 
any y > 1/2, there is a second order phase transition (in 
the DP universality class) to a surviving active phase at 
z = z c (y). Finally, when y = 1/2 the cluster dies still 
with probability one, but the average life time is infinite 
when z = 1. This is then the case of compact DP [13J. 

In 2+1 dimensions we take a square lattice in space, 
and each site (x, t + 1) can be infected by any of the four 
sites (x ± ei, t), (x ± e2, t). Assuming still that at least 
one of these sites has to be active in order to activate 
(x, t + 1), we have now four wetting ("infection", "acti- 
vation") probabilities Pk, k — 1,2,3,4. Site percolation 
corresponds to Pk = p for all k, while bond percolation 
corresponds to Pi = p, Pk+i = P k + (1 — Pk)p for k > 1. 
These formulas are easily understood: In site percolation 
it does not matter how many of the neighbours are active, 
since the site will be wetted anyhow, if it can be wetted 



at all. In bond percolation, the chance to be wetted by 
k + 1 neighbours is equal to the chance that the first 
k of them succeeded, plus the chance that the last one 
succeeds if the first were unsuccessful. We have thus 4 
control parameters, but generically it will be sufficient to 
study a 2-dimensional subspace. We choose (somewhat 
arbitrarily) the following 2-parametric ansatz 

Pi = p, P k+1 = P k + (1 - P k )pq for k > 1 (5) 

with < p < 1 and q such that all Pk are positive and 
less than 1. This ansatz includes both site DP (q = 0) 
and bond DP (q — 1). When < q < 1/p, it has a similar 
interpretation to the one given for bond DP: if the first 
k attempts were all unsuccessful in wetting the site, the 
chance for the next one is not p but is qp. But Eq.JSJ) is 
independent of this interpretation and is legitimate also 
for q < 0, provided < Pk < 1 for all k. 

Simulations of this model are straightforward. At each 
time t we have two lists L id and L new of sites, together 
with an integer array S of size L x L, where L is the linear 
size of the system. At the beginning we put t = 0, S and 
L nGW are empty, and I/ id contains the coordinates of the 
seed (we use helical boundary conditions, i.e. each site 
is labelled by a single integer i with i = i + L 2 , and its 
four neighbours are i ± 1 and i± L. When going from t 
to t + 1, we first go through the list L \^ and increase S[i] 
by one unit, whenever i is a neighbour of a site in L id- 
At the same time we write each of these i into the list 
L ncw . In a second pass, we first erase L id) then we check 
which sites in L now are actually wetted (for this, we use 
the information stored in S). Those who are wetted are 
written into L Q id, and S[i] is reset to zero for all checked 
sites. We then erase L ncw , and are ready to go from t + 1 
tot + 2. 

There are two non-trivial tricks for improving this al- 
gorithm. The first is needed in the first order transi- 
tion regime. There, the survival probabilities decay ex- 
tremely fast with time t. Therefore we use enrichment, 
implemented recursively as in the PERM algorithm . 
Essentially, this makes a copy of every cluster which sur- 
vives up to a t where the fraction of surviving clusters is 
below some threshold, and lets the copy and the original 
evolve independently. 

The second trick is related to histogram reweighting 
pol l2l| , but it is done on the fly as in . Assume a 
cluster contains rik sites wetted by k neighbours (fc = 
1, ...4) and bk boundary sites which were not wetted, 
although they had k wet neighbours. Such a cluster is 
included in the sample with probability 

4 

P({n k }, {b k };p, q) = J] K k (1 - P*) bk - (6) 

k=l 

For any observable A, the average is 

(A) p , q - A({n k },{b k })P({n k },{b k };p,q) 

{n k },{bk} 

= M" 1 £ Adukjdh}), (7) 

events 
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where the first sum runs over all possible configurations 
and the second sum runs over all M actually simulated 
clusters. The average for some other pair (p' , q') of con- 
trol parameters is then given by 



(A> p /,,/ = M- 1 Y, A({n k },{b k }) 



P({n k },{b k };p',q') 
P({n k },{b k };p,q) ' 



The ratio of probabilities is actually calculated by mul- 
tiplying the corresponding ratios of P k resp. 1 — P k for 
each site which is wetted (resp. not wetted) during the 
build-up of the cluster. 

This is useful, if we want to estimate several averages 
during the same run. It is particularly useful, if we want 
to estimate derivatives of (A) p q with respect to p or q. 
The latter is needed for checking scaling laws, as we will 
discuss below. Such derivatives can be obtained directly 
from Eq.@, e.g. 



with 



~fil{A)p,q = M~ x Y A({n k },{b k })W 

events 



w = J2 ojWM 



(9) 



wetted sites 



J2 QjMi-PkM)- (io) 



non— wetted sites 



We made several types of runs. Most extensively, we 
started from a single seed and measured the averages of 
the number N(t) of sites wetted at time t, of the sur- 
vival probability P(t), and of the rms. distance R 2 (t) of 
wetted sites from the seed [23j. In several of these runs 
we also calculated reweighted averages of N(t) and/or 
d(N(t))/dp. In addition, we made also some simulations 
with fully active initial conditions, where we measured 
the decay of the density pit) of sites active at time t. 
Finally, we made some runs also with initial conditions 
where half of the lattice was fully active and the other half 
was empty [T(ij| . In that case one has two phase bound- 
aries, and one can watch how these boundaries evolve 
with time. 

In all cases, lattice sizes were sufficiently large that we 
have no finite size corrections. This is most easily checked 
in single-seed runs: One simply has to check that the 
wetted region never hits the boundary of the lattice. In 
the other cases the check is less straightforward, but we 
always verified that L/2 3> t x l z , where z is the dynamical 
exponent. In most cases this implied L > 1000 for t w 
40, 000. 



III. THEORETICAL PREDICTIONS 

For each value of g, we expect that there is a critical 
value p c (q) such that the activity dies out for p < p c {q), 
while lim^oo P(t) > for p > p c (q). For all values of q, 
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FIG. 1: Phase diagram for the generalized DK model in 2+1 
dimensions. Below the transition line, activity dies out, while 
an active phase exists above it. The transition is in the DP 
universality class to the left of the tricritical point marked 
"TP", while it is first order to its right. The actual locations 
of the curve and of the tricritical point were obtained from 
numerical simulations. 



p c (q) is a decreasing function. Known numerical values 
are for q = and q = 1: p c (0) = 0.34457(1 Wsite DP 
0) and p c (l) = 0.2873381(1) (bond DP JUIH). The 
phase diagram should thus look roughly as in Fig. ^ At 
the tricritical point TP: (p,q) — (p*,q*), the observables 
N(t), P(t), R 2 (t), and p(t) should follow power laws 



N(t) ~tv(l + 0(r A )), P(t)~ 
R 2 (t)~t 2 /*(l + 0(t- A )), p(t) 



t- 6 (1 + 
-t~ tf (l + 



o(r A )), 
or A ))(n) 



Here, A is the exponent of the leading correction to scal- 
ing term. 

The upper critical (spatial) dimension is d c = 3. 
For d > d c one has the mean field exponents r\ = 
Q,S / = l,z = 2, and 5 = 1/2. The predictions of the 
e— expansion with e = 3 — d are 



Tj = -0.0193e+O(e 2 ) p 
£' = l-0.0193e + O(e 2 ) 
z = 2 + 0.0086e + O(e 2 ) 



6 = - - 0.468e 
2 



0(e 2 



-0.019 
w 0.981 
« 2.009 

« 0.032 



(12) 



Here the numerical values are those obtained for d = 2 
by simply neglecting all terms higher than linear in e. 
These values satisfy, up to terms 0(e 2 ), the hyperscaling 
relation [TlTl2l| 



n 



+ 5 + 5' = d/z. 



(13) 



To derive the latter, remark that P(t) is the probability 
that there exists at least one path of active sites con- 
necting the site (x , t ) — (0, 0) to any of the sites (x, t). 
Similarly, p(t) is the chance that there exists at least one 
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path connecting any site (x, 0) to (Q,t). The probabil- 
ity that (0, 0) is connected to (0, t) scales, according to 
our Ansatze for N(t) and R 2 (t), as P conn - N(t)/R d ~ 
fn-d/z _ Q n ^ ne th er hand, if there is basically at most 
one percolating cluster near x = (an assumption which 
breaks down for d > d c ), and if the site (x, 0) belongs 
to it, then the condition for (0,t) to be wetted is the 
same as that for (0, t) to be connected to (0,0). Thus 
-Pconn ~ P(t)p(t), which gives immediately Eq. l(T5)l. 

In the vicinity of the tricritical point one should dis- 
tinguish between the behaviour on the critical curve (but 
with q q*) from the behaviour off the critical curve. 
For the former, we expect for each observable a scaling 
ansatz with scaling variable (q* — q)t x , e.g. for N(t) we 
have 

N(t;q,p=p c (q))=t*>F((q*-q)t x ). (14) 

Notice that it would be more standard to use {q*-q) 1 / x t 
as scaling variable, but then the scaling function replac- 
ing F(z) would be singular at the origin. The advantage 
of Eq. I|14f) is that F(z) is analytic at z — 0. On the other 
hand, for q = q* and p ^ p c {q*) we have |27| 

N(t;q*,p)=t*>G((p c (q*)-p)t y ). (15) 

Again, the standard way of writing this would be in terms 
of the scaling variable (p c (q*) — p)" ll t with v\\ = 1/y, 
the advantage of Eq. I|15|> being that G(z) is analytic at 
z = 0. 

Eq. l|15f) shows that the correlation time scales as r oc 
(Pc(q*) — p)" 11 > and therefore the correlation length scales 
as £ oc (p c {q*) — p) u± with v±_ = v\\/z. The cross over 
exponent <p is finally defined as 

<t> = x/y. (16) 

The e— expansion gives 0,13 

= l + 0.0193e+O(e 2 ) « 1.019, 

cj) = --0.0121cm 0.488, 
2 

z/i = - +0.0075e w 0.507. (17) 

Since the tricritical point TP is a fixed point of the 
renormalization group flow which is instable in both di- 
rections, we should slowly cross over to the scaling laws 
for normal DP, if we are on the transition line with 
q < q* . Thus, normal DP asymptotics is expected on 
the entire transition line left of TP. For DP one has the 
same scaling laws Eci. ()ll|l . but with different exponents 
IM IH: V = 0.2303(4), 8 = 5' = 0.4509(5), and 
z = 1.7666(10). These values satisfy the same hyperscal- 
ing relation Eq. (|13fl . 

The identity 8 = 8' follows from a special time reversal 
symmetry which holds for bond and site DP, but not for 
the DK model in general. As we said above, both P{t) 
and p(t) are probabilities that one point at the boundary 
of a time slice < t' < t is connected to some point on 



the opposite boundary. For bond DP one sees easily that 
both are the same, i.e. p(t) = P(t) exactly. For site DP 
one gets the same, if one assumes also all sites with t' = 
and t' = t as wettable with probability p, which gives e.g. 
P(l) = p(l) = p[l - (1 -p) 4 ]. No relation like that holds 
for the general DK model. 

On the part of the transition line with q > q* one has 
compact clusters with only few small holes and a slowly 
moving rough surface. This is analogous to critical com- 
pact DP, except that the latter occurs (in 1+1 dimen- 
sions) only at a single point and has clusters without any 
holes at all. Mean field arguments ^} suggest that the 
transition is first order along this part, in the sense that 
the density of active sites in the stationary state is dis- 
continuous. The growth of a cluster is then similar to the 
growth of a droplet in a usual first order phase transition, 
except that now one of the two phases is not fluctuating. 

Consider now an initially straight interface between the 
two phases. If p < p c {q), this interface will move into the 
direction of the active phase, i.e. the inactive phase will 
win. The opposite is true for p > p c (q), i.e. the value 
°f Pc(q) is fixed by the condition that neither phase will 
win in the long time limit Notice, however, 

that this does not mean that the interface at p = p c (q) 
will not move during short times. Due to the asymmetry 
between the two phases, we expect the generic behaviour 
of non-equilibrium ("growing") interfaces to apply, i.e. 
the velocity of the critical interface should scale with time 
as v ^ i -2 / 3 for large t |1S||. 

Let us finally discuss the growth of a cluster start- 
ing from a single site for q > q* ,p — p c (q). As in the 
Grassberger-Chate-Rousseau model with negative par- 
tial immunization, it is easier to survive for the clus- 
ter in its already occupied region than to spread out 
from this region. Although the details of the two mod- 
els are different, we might thus expect that we find 
also in the present model a stretched exponential law, 
P(t) ~ exp(— ct a ), for basically the same qualitative rea- 
sons as in the Grassberger-Chate-Rousseau model [Trjj . 
We have however no detailed prediction on the exponent 
a. 



IV. RESULTS 

A. Gross features of the phase diagram 

For ordinary DP, the observable most sensitive to the 
precise value of p c is N(t). We expect therefore that N(t) 
is also the best observable to locate the tricritical point. 
In Fig. [5] we show N(t) for twelve pairs of control param- 
eters (p,q). These represent three values of p (critical, 
sub-, and supercritical) for each of four values of q. The 
latter are chosen to be: (i) site DP (q = 0); (ii) in the 
crossover region from TDP to DP (q = 2.4); (iii) close 
to TDP (q « q* w 2.8) and (iv) q > q* . First of all, we 
see from Fig. [5] that there are indeed two straight lines. 
These are of course the candidates for ordinary and for 
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TABLE I: Estimates of p c (q). 
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FIG. 2: Log-log plots of N(t) for four different values of q, 
and for three values of p for each q. The three p— values are 
chosen to be subcritical, critical, and supercritical (bottom to 
top). Within each triple, they differ by less than 0.05%. The 
four q— values are DP, crossover from TDP to DP, TDP, and 
first order (top to bottom). 
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FIG. 3: Log-log plots of N(t) for several values of q, each at 
the best estimate of p c (q)- 



tricritical DP. Secondly, we see that we can determine 
p c (q), for each value of q, with rather high precision. For 
q = q* we just have to look for a power law, just as for 
q <C q* ■ For q in the cross-over region, i.e. close to q* 
but not at q* , the determination of p c (q) is less easy. But 
we still can get rather precise estimates for q slightly be- 
low q* , if we assume that the slopes of the critical curves 
approach the slope r\ = 0.2303 of the critical DP curve 
monotonously. For 5> q* , finally, it is more easy to de- 
termine p c {q) from the requirement that the velocity of a 
straight interface scales as i~ 2 / 3 (see details below) than 
from the behaviour of N(t). Numerical results for p c (q) 
obtained in this way are collected in Table 1. Several 
critical curves are plotted again in Fig. [3] 
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FIG. 4: Log-linear plots of t v N(t) for several values of q i 
and p ~ p c (q), with 77 = —0.353. 



B. Scaling at the tricritical point 

Plots like those in Fig. [2] are of course not suitable 
for precise determinations of critical parameters. For a 
closer look at the tricritical region, we show in Fig. 0] 
the product t"^N{t) with r] = -0.353. This value of 
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i] is our best estimate, i.e. the tricritical curve should 
be horizontal if there were no corrections to the scaling 
behaviour. Actually we see that there are quite strong 
corrections, manifesting themselves as a bump at t ~ 3. 
The data shown in Fig. [2] and similar data at other near- 
by values of q show that 

q* = 2.792(6), V c = 0.1831534(5) + 0.458(2.792 - <f ), 

(18) 

and 

X] = -0.353(9). (19) 

Although the last estimate has the same sign as the pre- 
diction from the e— expansion, it is larger by nearly a 
factor 20! 

Once we have fixed the tricritical point, we can imme- 
diately obtain 5' and z from the scaling of P(t) and R 2 (t). 
We do not show the corresponding data. The estimates 
found by standard extrapolation methods are 

8' = 1.218(7), z = 2.110(6). (20) 

For z we find that the deviation from the mean field 
value z — 2 has the same sign as predicted (TDP spreads 
subdiffusively, in agreement with the intuitive notion that 
spreading should be slowed down compared to ordinary 
DP). But again the difference from mean field is an order 
of magnitude larger than predicted. For the correction 
to 8' , the e— expansion predicted even a wrong sign. 

Finally, in order to measure the exponent 8, we made 
runs with fully active initial conditions. Again we do 
not show the raw data, as they contained little surprises. 
They give 

5 = 0.087(3). (21) 

In spite of the huge difference from mean field behaviour 
(where 8 = 1/2), this is in surprising good agreement 
with Eq. i|12|) . The hyperscaling relation is very well 
satisfied with these estimates, 

r? + 8 + 5'-d/z = 0.004(13), (22) 

showing that they are at least internally consistent. 

C. Scaling near the tricritical point 

A popular method to find correlation length and cross- 
over exponents are data collapse plots. In view of a scal- 
ing law like Eq. (|15fl . it seems natural to plot N(t)/t v 
against (p c (q*) — p)t v and determine y (or, equivalcntly, 
z^ll) so that all data fall onto a single curve. Such a plot is 
shown in Fig. EI where we have chosen y = 0.87. It seems 
to give a perfect collapse, i.e. there are no indications of 
any scaling violations. This might seem strange in view of 
the rather strong violations seen in Fig.0] It is true that 
we suppressed them by plotting only data with t > 30, 
but this should not eliminate them completely. Rather, 
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FIG. 5: Log-linear plots oit~ v N(t) against (p c (q*)-p)t y , with 
q* = 2.8, t) = 0.362, and y = 0.87. Only points with t > 30 are 
shown. Values of p range from 0.17693 to 0.185350. Notice 
the seemingly perfect data collapse in spite of the slightly 
wrong tricritical parameters and the scaling violations seen 
in Fig. H 



they are not seen in Fig. \5\ because of the coarse scale 
on the y-axis, and because we also changed slightly the 
tricritical parameters. We used q* — 2.8 and r\ = 0.362, 
while the plot would be definitely less clean if our best 
estimates had been used. 

A much more reliable and precise method consists in 
measuring dN(t; q,p)/dp. As we have pointed out in 
Sec. 2, this can be done straightforwardly. From the scal- 
ing ansatz Eq. (|15f) we have 

aiog^(y, Pe(g ')) =fy diogG(z) 

op dz 

up to corrections to scaling. This quantity was found to 
depend very weakly on the precise value of p. Thus the 
main uncertainty comes from its dependence on q, to- 
gether with our rather large error for q* . We plot there- 
fore in Fig. El the quantity t~ v d\ogN(t)/dp for various 
values of q close to q*, taking for each curve our best 
estimate of p c . We see quite strong corrections to scal- 
ing (as we expected), but they do not prevent us from 
a very precise estimation of y. Fitting with a correc- 
tion to scaling exponent A = 0.66, i.e. with an ansatz 
d\ogN(t)/dp x t y (l + const/i 66 ), we obtain 

y = 0.865(3), v\\ = 1.156(4), (24) 

and, using the previously determined value of z, v±_ = 
0.547(3). Using this, we obtain also = v\\6 = 0.100(4) 
and f3' = Un8' = 1.408(10), where f3 and j3' are defined via 
the scalings lim t _> 00 p(t) ~ (p - p c {q*)) p , lim t ^oo P(t) — 
(p " Pc{q*)Y j both for p > p c . Again these estimates, 
albeit being close to their mean field values, deviate from 
them much more than predicted by the e— expansions 
(which is particularly simple for which should be equal 
tol + (9(e 2 )). 
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FIG. 6: Log-linear plots of t v d log N(t)/dp against t. 
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FIG. 7: Log-linear plots of t~ v N(t; q,p c (q)) against (q* —q)t x , 
with q* — 2.792, r\ = 0.353, and x = 0.31. Only points with 
t > 24 are shown. Values of q range from 1.6 to 3.5. 



Finally, in order to measure the cross-over exponent 
cj>, we can either use again a data collapse plot, or we 
can try to estimate the derivative of N(t; q,p) (or of any 
other observable) along the curve p — p c (q)- The latter 
seems however not easy. It is of course straightforward 
to measure the derivative in any given direction. But 
the orientation of the transition line p — p c (q) is only 
approximately known, and the derivative should depend 
strongly on this orientation. Thus we determined 4> only 
from its data collapse, in spite of the shortcomings men- 
tioned above. 

Indeed, the situation is now much worse than in the 
data collapse plot for , partly because of the un- 
certainties in p c (q)- In Fig- El we show the values of 
t~ v N(t;q,p c (q)) plotted against (q — q*)t x . The values 
of 7/ and of p c (q) are as determined above, and x = 0.31. 
The data collapse is now much worse than in Fig. [S] This 
may be due to the fact that Eq. I|14|) has even larger cor- 



rections to scaling than Eq. (|15(l . but it may also be due 
to errors in estimating p c {q). The data collapse would 
improve considerably, if we would shift p c (q) for q < q* 
systematically towards higher values. In that case the 
critical curves In N(t) versus lnt would no longer be con- 
cave, i.e. din N(i)/d\nt would no longer approach ap- 
proach the value t^dp = 0.2303 monotonically from be- 
low. Although we cannot rule out such a behaviour, we 
prefer to keep the estimates of p c (q) and to increase the 
estimated error of cf>. We thus obtain x = 0.31(3) and 



= 0.36(4). 



(25) 



Again the deviation from the mean field value (cj) = 1/2) 
is in the same direction as predicted by the e— expansion, 
but is much larger in absolute value. 

The numerical values of the tricritical exponents, to- 
gether with the previous estimates from the e— expansion 
and from the simulations in Ref . fT(il | , are collected in Ta- 
ble 2. 



D. The first order transition 

As we said, data like those shown in Figs. 12131 are 
strongly suggestive of stretched exponentials, but it is 
notoriously difficult to obtain the precise asymptotic be- 
haviour from such curves. It seems not even clear which 
of the curves for q = 3.2 in Fig. 121 is closest to the criti- 
cal one. An alternative way to obtain the critical point 
in this case consists in the following. We start with ini- 
tial conditions where one half of the square lattice of size 
L x L (let us say all sites with L/4 < x < 3L/4) is 
active and the other half is dead. Thus we have two in- 
terfaces, at x — L/4 and at x — 3L/4 (we assume that 
L is a multiple of 4, and boundary conditions are peri- 
odic in the x direction). We then measure the density 
p(x,t) as the system evolves. The initial step functions 
will become smooth sigmoidals, corresponding to fuzzy 
interfaces. Their positions are measured by measuring 



X(t) 



(L-l)/2| 



L/4. (26) 



Initially, X(t = 0) = 0. For large t, X(t) increases lin- 
early with t if p > p c (q), while it decreases if p < p c {q)- 

Data for q = 4.0 are shown in Fig. [H] where we have 
plotted lnA(t) against \nt for four different values of 
p. In addition we show in this figure the power law 
X(t) ot i 1 / 3 which is generically expected for a rough 
non-equilibrium surface which neither grows nor recedes 
for t — > oo. We see that this is indeed the asymptotic 
behaviour for p — 0.137830(3), which is thus our best 
estimate for p c (4.0). 

Similar behaviour is found also for other values of q > 
q* , except when q is very close to q* . There the transient 
behaviour seen in Fig. [S] for t < 20 extends to much 
larger t, reflecting the increased internal fuzziness of the 
interface. 
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TABLE II: Estimated critical exponents for tricritical DP in 2+1 dimensions. Values for the spatial fractal dimension Df 
and for the exponents 7, a, and r/± not discussed in the text were obtained by means of scaling an d hy perscaling relations as 
indicated in column #1. They are included here for easier comparison with the simulations of Ref. |l0|) (last column). 
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FIG. 8: Log-log plots of X(t) versus t for q = 4.0 and four 
values of p. The straight line corresponds to X(t) oc t 1 ^ 3 , the 
generic behaviour of a non-equilibrium surface which neither 
grows nor recedes asymptotically. 
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FIG. 9: Log-linear plots of R 2 (t)/t - 7 versus t, for q = 4.0 and 
four values of p. 



It is an interesting question how clusters starting from 
a single site grow at the first order transition line. A 
priori one might expect the behaviour to be similar to 
that of a small droplet in a gas at boiling temperature, 
or a small domain of inverted spins in an Ising magnet at 
T < T c . But these analogies might break down due to the 
essentially non-equilibrium nature of DP. Let us concen- 
trate again on q = 4.0 (at other values of q, the behaviour 
was similar). Quantities R 2 {t) and N{t)/P{t) (i.e., the 
average size of surviving clusters) are plotted against t 
in Figs. and Actually, in order to see the scaling 
more clearly, in both cases the plotted quantity was first 
divided by a suitable power of t, so that the curves would 



be straight horizontal lines if there were pure power laws. 
But the curves, in particular those for the critical value 
p = 0.13783, are far from straight. They are compatible 
with approximate scaling laws R 2 ~ N(t)/P(t) ~ t - 5 
to i ' 6 , but we are still far from the asymptotic regime. 
For t — > 00 we should expect that surviving clusters are 
roughly circular and have a compact interior with den- 
sity p = N(t)/P(t)/R d (t). If the transition is indeed first 
order, this density should tend to a positive constant for 
t — > 00. As seen from Figs. I§1 and ITU1 this is indeed the 
case, although this limit is reached rather late. 
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FIG. 10: Log-linear plots of N(t)/ P(t)/t - 7 versus t, for q = 
4.0 and four values of p. 



V. DISCUSSION 

The main result of the present analysis is the rather 
poor agreement with the e-expansion. This is somewhat 
surprising, since the upper critical dimension is d c = 3 
for this problem, and therefore terms linear in e should 
give a rather good description in two dimensions. It is 
even more strange, since the contributions of order e are 
extremely small for all critical exponents - except for <5, 
and it is only for 8 that there is fairly good agreement. 
At least, for most exponents (except 5') the deviations 
from the mean field predictions have the same sign as 
predicted to order e. 

Our results are also in rather poor agreement with 
those of the simulations of Liibeck E.g., he found 

4> = 0.55(3) which is 5a from our estimate, and which 
deviates from the mean field value 1/2 in the opposite 
direction than the O(e) term. His estimate v± = 0.59(8) 
agrees with ours within the error bars, but his estimate 
f3 = 0.14(2) exceeds ours by 2 standard deviations. The 
reason for this discrepancy is not clear. It could be due to 
the fact that the analysis in [lfj was mainly based on data 
collapse plots, which makes it very difficult to take into 
account possible corrections to scaling. As we pointed 
out in Sec.4.C, such plots can hide even large systematic 
deviations from the supposed scaling behaviour. 

Overall, the simulations presented in this paper pro- 
duced hardly any difficulties or surprises. Except for the 
data collapse for the cross-over exponent <j>, all correc- 
tions to scaling were rather modest. We thus believe 
that our analysis is robust and does not hide any large 
systematic errors. 

This is in marked contrast to the PCPD [2|, which can 
be viewed as an alternative tricritical version of stan- 
dard DP. It would be nice to have a model where the 
PCPD arises as a singular limit of TDP. The version of 
the generalized Domany-Kinzel model which we used in 
the present paper does not give rise to the PCPD in any 



limit. 

After having discussed tricritical directed percolation, 
we should mention also tricritical ordinary (undirected) 
percolation. A field theoretic renormalization treatment 
of this problem has been given recently by Janssen and 
coworkers |(| . Discrete lattice model versions can be for- 
mulated in close analogy to the present work. Let us 
consider the dynamic version, i.e. the generalized epi- 
demic process (GEP) |33|. This is very similar to directed 
percolation viewed as an epidemic process, except that 
lattice sites which had once been infected cannot be re- 
infected again (this is also known as the SIR model, for 
"susceptible-infected-removed"). In a generalized GEP 
(GGEP) one can modify the probabilities Pk to be in- 
fected by k infectious neighbours in the same way as 
we did in the present paper (Eq. 0, but one can also 
use other Ansatze for Pk with the same qualitative be- 
haviour. Typically, a tricritical regime is reached when 
Pk increases sufficiently fast with k. Models of this type 
were studied some time ago by Cieplak, Robbins, Koiller 
and others [33] • A comparison of their results with the 
theoretical work of |(J would be very welcome. Also, 
there are some surprising claims in |32|. e.g. it is claimed 
in some of these papers that there is also a transition be- 
tween pinned fractal (i.e. percolation like) clusters and 
compact clusters with self-affine surfaces in 1+1 dimen- 
sions. One should expect such a transition to occur in 
higher dimensions, but like the transition to compact DP 
in the 2-d Domany-Kinzel model, it should occur in 1+1 
dimensions only at the boundary of the control parame- 
ter space. 

Finally, it is instructive to compare different epidemic 
models, and we shall finish this paper with a short review 
of how the different schematic models discussed in the lit- 
erature are related. Let us first assume we can have three 
different types of individua in a population: Healthy sus- 
ceptible ones (S), ill and infective ones (X), and "re- 
moved" individua (R), which might either mean immune 
or dead. Then we have the following basic schemes: 

• The only reaction is S + X — > X + X, i.e. ill in- 
dividua remain infective forever. This gives Eden 
growth |Tsj . 

• In addition we have X — ► R, i.e. ill individua have 
a finite life time and are removed thereafter. This 
gives the GEP. 

• Instead of X — > R we have X — > S, i.e. after indi- 
vidua have recovered they become again suscepti- 
ble. This gives DP. 

• If we add a reaction R — > S to the GEP, it should 
bring the model into the DP universality class, ex- 
cept when the rate of this last reaction is very small. 
In the limit when this rate is much smaller than all 
other rates, we obtain the Bak-Chen-Tang (3^ for- 
est fire model. 
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If we have in addition a fourth class (W) of individua 
weakened by the contact with an infective, but neither ill 
nor infective themselves, then we get three more schemes: 

• If we change the recovery in DP to X — > W and 
add W + X — > X + X , i.e. recovered individua have 
a different chance to be re-infected by a subsequent 
contact with ill ones, we obtain the model of Grass- 
berger, Chate, and Rousseau 0|. Notice that in 
this model weakened individua stay weak forever, 
and never regain their strength. 

• If we include in addition the reaction W — > S, then 
weakening is only transient and we should expect 
to be again in DP. 

• If, however, transient weakening happens not after 
but instead of an infection, i.e. if we add S + X — + 
W + X, W + X -» X + X, and W -> S to DP, we 
obtain the model discussed in the present paper. 

• If we add S + X -> W + X and W + X -» X + X 



(but not W — ► S) to the GEP, we end up at the 
GGEP. 

Obviously this list is not exhaustive. It seems however 
to contain all interesting cases, if we assume that the 
reservoir of susceptible individua is unlimited (conserva- 
tion of the total population size is irrelevant), and that 
there is some local saturation mechanism which prevents 
any of the densities to explode. Whether this is indeed 
true, or whether there exist more such models with in- 
teresting novel behaviour, is still an open question. 
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